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Abstract. Simple analytical parametrizations for the ground-state energy of the one- 
dimensional repulsive Hubbard model are developed. The charge-dependence of the 
energy is parametrized using exact results extracted from the Bethe-Ansatz. The 
resulting parametrization is shown to be in better agreement with highly precise data 
obtained from fully numerical solution of the Bethe-Ansatz equations than previous 
expressions [Lima et al., Phys. Rev. Lett. 90, 146402 (2003)]. Unlike these earlier 
proposals, the present parametrization correctly predicts a positive Mott gap at half 
filling for any U > 0. The construction is extended to spin-dependent phenomena by 
parametrizing the magnetization-dependence of the ground-state energy using further 
exact results and numerical benchmarking. Lastly, the parametrizations developed 
for the spatially uniform model are extended by means of a simple local-density-type 
approximation to spatially inhomogeneous models, e.g., in the presence of impurities, 
external fields or trapping potentials. Results are shown to be in excellent agreement 
with independent many-body calculations, at a fraction of the computational cost. 
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1. Introduction 

The archetypical strong-correlation phenomenon is the Mott insulator [Ij. It is well 
known that the first-principles description of such systems by means of density- 
functional theory (DFT) encounters severe difficulties. The single-particle (Kohn-Sham) 
gap calculated with standard DFT methodology is not the same as the many-body gap, 
even in principle and if no approximations at all are made during the calculation. In the 
particular case of the Mott insulator, it is known that a proper description of the Mott 
gap is obtained by adding to the single-particle Kohn-Sham gap a correction arising 
from the derivative discontinuity of the exchange-correlation functional |2H7J. 

Many modern density functionals do have an implicit derivative discontinuity due 
to their orbital dependence, but this affords at best an incomplete description of the 
Mott state and the Mott gap [3],[6l[7] . Only very few density functionals have an exphcit 
discontinuity as a function of the density, among them the 2-electron functional of 
Mori-Sanchez, Cohen and Yang [3] which is based on the so-called flat-plane condition, 
devised by the same authors, [3H5] and the Bethe-Ansatz LDA (BALDA) of Lima et 
al. [HIIH], which is based on an approximate analytical parametrization of the Bethe- 
Ansatz solution for the ground-state energy of the one-dimensional Hubbard model 
that becomes exact in several important limits. Both of these functionals do properly 
account for the Mott insulator. Only the BALDA, however, has been parametrized in a 
way that allows its application to a wide range of physical systems. As a consequence, 
BALDA and variations thereof has been applied to inhomogeneous correlated many- 
electron systems as well as to correlated many-atom systems in optical lattices and 
traps [9HT8]. 

However, in these applications it has become clear that the parametrization by Lima 
et al., also referred to as LSOC parametrization, after the initials of its developers, 
still does not provide a fully correct description of the Mott gap. In particular, at 
small U, the LSOC expression for has a derivative discontinuity, but the resulting 
Mott gap is negative, effectively predicting the Mott state to be a (strange) metal 
instead of an insulator. It is not trivial to correct this behaviour, as any change to the 
LSOC expression must preserve the exact limits and properties already built into it. 
Thus, instead of merely algebraic adjustments, it becomes necessary to understand and 
improve the physics missing from the LSOC parametrization in this regime. 

A key aspect of the Mott insulator is that it is a nonmagnetic state of matter, 
i.e., the insulating nature is not the result of antiferromagnetism. Therefore, the 
solution to the problem just described must be expressed, within DFT, in terms of the 
charge- density only, and cannot make use of spin densities and spin-density-functional 
theory (SDFT). On the other hand, many correlated systems do have magnetic phases. 
Therefore, a more complete description of strong correlations must properly account for 
both, the insulating state and various types of magnetic states, as well as their possible 
coexistence. While such description is possible within Bethe-Ansatz based SDFT by 
performing a fully numerical solution of the spin-resolved Bethe-Ansatz equations and 
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interpolating between the resulting data points every time the exchange-correlation 
energy needs to be evaluated, this procedure is very inconvenient for Kohn-Sham 
calculations, where the functional must be evaluated hundreds or thousands of times 
during the iterations towards selfconsistency. A simple parametrization of the spin 
dependence would allow straightforward application of BA-DFT methodology to spin- 
dependent phenomena in electronic systems and hyperfine-label-dependent phenomena 
in optical lattices. 

The present paper reports progress along both of these lines. In Section[2]we identify 
a shortcoming of the standard parametrization used in BALDA and propose a simple, ad 
hoc but physically motivated, variation of it that is more accurate and whose Mott gap 
properly remains positive for all values of the interaction U. The revised parametrization 
also considerably improves the description of the metallic (e.g., Luttinger liquid) phases. 
In Section [3] we employ exact analytical results extracted from the Bethe Ansatz to 
further generalize this revised parametrization to spin-dependent situations. In Section 
m we use a simple local-density type approximation to extend our results to spatially 
nonuniform systems. Density-Matrix Renormalization Group (DMRG) and Lanczos 
calculations are performed to test and validate this approximation. Section |5] contains 
a brief summary. 

2. Improved description of the Mott gap 

The task at hand is to obtain a simple and accurate analytical approximation for the 
ground-state energy Eq of the inhomogeneous one- dimensional Hubbard model (IDHM), 



where rijo- = cl^Cicr is the spin-resolved particle-density operator at site i, and Ci„ are 
fermionic creation and annihilation operators, t is the hopping parameter (taken to be 
the unit of energy), U the on-site interaction and Vi an on-site potential that makes the 
system spatially inhomogeneous. In applications to electrons in crystal lattices, Vi can 
describe inequivalent atoms in the lattice, while in applications to cold atoms in optical 
lattices it accounts for the trapping potential. 

The commonly used expression for the per-site ground-state energy (eo) of the 
IDHM is the LSOC parametrization, given by [8] 



where n = n^- + n^. is the charge density and the interaction U enters eo(n, U) through 
the interaction function f3{U), which is determined from 



By construction, this expression becomes exact for ?7 — )■ and any n (where /3 = 2), for 
f/ — )■ oo and any n (where /3 = 1), and for n = 1 and any U (where < /3 < 1), and 
provides a reasonable approximation to the full Bethe-Ansatz solution inbetween. 




(1) 




(2) 




(3) 
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In the LSOC parametrization, the interaction function I3{U) is independent of the 
particle density and of spin. This independence is algebraically very convenient, as it 
allows one to determine (3{U) outside the selfconsistency cycle of DFT instead of having 
to recalculate it any time the charge (or spin) density changes. However, it is physically 
incorrect, as the relation between the bare interaction parameter U and the correlation 
energy must depend on the charge density, e.g. through screening. As a consequence 
of the density-independence of (3, LSOC retains the sinusoidal density- dependence of 
the energy, which is correct only aX U = and f/ — ?■ oo, also for all intermediate values 
of U. 

Thus, it becomes necessary to modify the interaction function in a way that allows 
it to change with the density. The additional density dependence, however must not 
spoil the behaviour in the three limits exactly recovered by LSOC. This condition severely 
restricts the modifications that could possibly be made to the LSOC parametrization. 
The form we adopt is 



where /3{n,U) = /3(f/)"*-"'^^ and a{n,U) = n^^^. We stress that this particular 
form has not been derived from first principles but is a physically motivated ad hoc 
modification designed to restore the density-dependence of the interaction function 
through the replacement f3{U) f3{n, U), while preserving all exact limits obeyed by the 
LSOC expression. The specific form chosen for the exponent a{n, U) is a consequence 
of the later generalization to spin-dependent phenomena, as explained in Sec. [3l 

Figure [1] compares the present parametrization (jlj) to data obtained from a fully 
numerical (FN) solution of the Bethe-Ansatz integral equations. For comparison 
purposes, the earlier LSOC parametrization (jj]) is also included. To distinguish the 
present from the LSOC expression, we label the curves corresponding to the former by 
FVC. As the insets show, the relative deviation of FVC data from FN data is typically 
less than 2%, and at most ~ 4%. This is the same size of error of the local-density 
approximation itself (quantified by comparing BALDA/FN data for small Hubbard 
chains to results from exact diagonalization) , so that to within the accuracy of the LDA 
the present parametrization is a faithful representation of the full BA solution for the 
entire parameter range, including values of f/ ^ 6t that cannot be realized in solids but 
occur in systems of trapped cold atoms. 

Our main motivation for developing Eq. (jlj), however, was the incorrectly negative 
Mott gap predicted by the LSOC expression between [7 = and U = 2t. The Mott gap 
Egap can be evaluated analytically from the expression for eo(n, U), either by explicitly 
calculating the derivative discontinuity or by taking total-energy differences [9]. From 
the LSOC parametrization one obtains [9] 




(4) 




(5) 
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Figure 1. Per-site ground-state energy as a function of: a) density, for U — 6t; b) 
interaction, for n ^ 0.5. Insets: percentage deviation from the fully numerical (FN) 
data, defined as 100(e^^^(^^o^) - e^^)/e^^. 



whereas the present parametrization leads to 
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As illustrated in Fig. [21 the additional terms arising from the present expression push 
the gap upward, avoiding that it becomes negative between [/ = and U = 2t as was, 
incorrectly, the case within the LSOC expression [H]. For U > 2t both the LSOC and 
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Figure 2. Mott gap Egap of the infinite homogeneous chain obtained from numerical 
solution of the full BA equations, compared to analytical predictions obtained from 
the LSOC and the present (FVC) parametrization. 

the present (FVC) gap are positive, the latter being significantly closer to the numerical 
BA results than the former, although neither reproduces the subtle nonperturbative 
behavior for f/ — )■ 0, in spite of the logarithmic term in Eq. (Q. Only within the latter, 
however, the Mott gap is everywhere positive. 

3. Extension to spin-dependent phenomena 

In spin-polarized situations, the energy depends on the spin density m = — ni, 
in addition to the charge density n = n-^ + and the interaction U. This 
dependence is not included in the LSOC parametrization, which can therefore not 
be applied to study spin-dependent phenomena for electrons in solids or hyperfine- 
polarization-dependent phenomena for atoms in optical lattices. We note that BALSDA 
calculations for the IDHM have already been performed in the context of cold atoms 
in optical lattices [T9H2T]. In lieu of an analytical parametrization, these works 
resorted to a fully numerical (FN) solution of the BA integral equations. An analytical 
parametrization would substantially simplify this approach, making the BALSDA as 
easily implementable as the BALDA one. 

An additional advantage of the analytical approach over the numerical one is that 
in the solution of the BA integral equations one cannot specify from the outset the 
density and magnetization of the system one is interested in. Rather, one has to specify 
the upper and lower limits of the integrals, and obtains the densities as part of the 
solution. This is clearly inconvenient for DFT, where the energies have to be evaluated 
as functions of the densities. More generally, it is desirable to be able to specify the 
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system under study in terms of physical observables, such as the densities n and m, 
instead of in terms of auxihary quantities, as the hmits of the BA integrals. 

Analytical expressions also permit one to derive further analytical results for other 
quantities. Our present derivation of closed expressions for the Mott gap is one example, 
and the analytical derivation and solution of Euler equations determining the phase 
diagram of harmonically trapped fermions on optical lattices |T7] is another. For all 
these reasons, a simple but reliable parametrization of eo{n,m,U) can be useful for 
various types of calculations, within SDFT and beyond. Therefore, we next present 
an analytical parametrization of eQ{n,m,U) and use it to construct a Bethe-Ansatz 
local-spin-density approximation (BALSDA) for the IDHM. 

In order to generalize the LSOC and FVC expressions to spin-dependent situations, 
we one more time follow the basic philosophy of constructing an analytical interpolation 
function recovering exactly known limits as function of the physically relevant variables. 
Several such exact results for eo{n, m, U) are known [221421] . For non- interacting systems 
{U = 0), 

4 /7m\ /vrmx 
eo(n, m, (7 = 0) = sin I — 1 cos I — — 1 . (7) 



71 \ 2 J \ 2 J 
For infinite interaction {U oo), 

2 

eo(n, m, [/ — )• oo) = sin(7rn). (8) 

vr 

For half- filled unpolarized systems {n = l,m = 0), 

Finally, for maximum magnetization {m = n), 

2 

eo{n,m = n,U) = sin(7rn). (10) 

71 

The LSOC and FVC parametrizations take m = 0, and express e(n, U) as a 
controlled analytical interpolation between ([H]), (E]) and the m = limit of ([7]). The 
construction of a more complete interpolation, recovering all four limits as functions of 
n, m and f/, is strongly constrained by these limits, but still not unique. We therefore 
impose five additional common-sense criteria: (i) avoid high-order polynomials, which 
can produce unphysical wiggles, (ii) avoid unusual special functions, and (iii) keep the 
form similar to the LSOC parametrization. This third condition is useful because by now 
the LSOC parametrization has been implemented and used by many groups [SHIl], so it 
will be easier to update to a new parametrization which has a similar form to the old one. 
However, as we have argued above, the LSOC expression does not properly describe the 
density-dependence of the ground-state energy at intermediate U . Therefore, we also 
require, as condition (iv), our spin-dependent generalization to reduce to the present 
Eq. (jl]) for m = instead of to Eq. ([2]). The final, fifth, additional condition is motivated 
by computational efficiency and is explained in Sec. |4]below. Even with these additional 
common-sense criteria, the form of the parametrization is not uniquely determined, and 
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a very large function space can still be explored. The particular choice made below 
was obtained by starting from the LSOC form and then building in, one-by-one, the 
additional exact limits and criteria. Many different variations have been explored, but 
we stopped when arriving at one whose deviation from the fully numerical solution 
of the Bethe-Ansatz equations was less than the typical error bar of the local-density 
approximation for this type of system. At this point, further improvements in the 
form of the parametrization become indistinguishable, in applications to inhomogeneous 
systems, from the intrinsic error of the LDA. 

All four exact limits and five supplementary conditions are incorporated by the 
choice 

2/3(n,m, U) 



n, m, U) 



sm 



where 



(3{n,m,U) = (3{U) 
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Here, P{U) is the same quantity employed in the LSOC parametrization. For zero 
magnetization, a{n,m,U) reduces to a{n,U) used in Eq. (jl]). Equation ( ITT]) is valid 
for ?7 > and n < 1, but it can be extended to n > 1 and to t/ < by standard 
particle-hole transformations [9|[T7 t l22 | l2l]. 
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Satisfaction of condition f llOp is of particular interest because this condition has 
no counterpart in the spin-independent situation. Moreover, it estabhshes a couphng 
between the charge and the spin dependence. Physically, maximum spin means that 
the Pauli principle keeps all fermions maximally apart for any U, just as infinitely 
repulsive interactions {U — oo) do for any degree of spin-polarization. For n = m the 
spin-dependent parametrization must thus reduce to the same limit as for f/ oo and 
recover the value /3 = 1. On the other hand, for m = the expression should reduce to 
the earlier form (jl]). This double requirement is the explanation of the particular form 
chosen for a{n,m, U) of Eq. ( !T^ and, consequently, for a(n, U) of Eq. (jl]). 

Figure [3] presents a comparison of the spin-dependence of this parametrization 
with data obtained from a fully numerical solution of the BA integral equations 
for intermediate parameter values, where the expression is not exact already by 
construction. Clearly, the spin-dependence of the homogeneous system is recovered 
to within the same precision as the charge dependence. 

4. Local approximation for inhomogeneous systems 

The main use of DFT and SDFT is in calculations for spatially inhomogeneous systems, 
where translational symmetry is broken and the density becomes position dependent. 
In the case of lattice models, inhomogeneity means that not all sites are equivalent. 
By means of the L(S)DA prescription, expressions for the energy of the homogeneous 
system, such as those described in the preceding section, can be used on a site-by-site 
basis to approximate the corresponding energy of the inhomogeneous system. 

Explicitly, the local-spin-density approximation to any energy component E is given 



where L is the number of lattice sites, and e = Wthl^ooE^L) j L is the per-site 
energy of the homogeneous system. This expression approximates the energy of the 
inhomogeneous system, where Ui and vary from site to site, by evaluating the energy 
density of the homogeneous system site by site at the densities of the inhomogeneous one. 
In DFT, including model DFT, this prescription is usually applied to the correlation 
energy, which for the Hubbard model can be defined as E^ = Eq — Emf, where 
Emf is the mean-field approximation to the ground-state energy Eq. Since Emf is 
simple to obtain, the task to approximate the correlation energy Ec[ni,mi,U] of the 
inhomogeneous Hubbard model, is thus reduced to that of approximating the per-site 
ground-state energy of the homogeneous one, eo(n, m, U). This is the quantity that we 
extracted above from the Bethe-Ansatz equations. 

The minimization of the resulting energy functional is conveniently carried out 
via selfconsistent Kohn-Sham (KS) calculations of the charge and spin densities 
and of the resulting ground-state energies. In such KS calculations the correlation 
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potentials (obtained by differentiating the correlation energy with respect to n and 
m, or, equivalently, to n-^ and n^) are evaluated once in each of the iterations 
of the selfconsistency cycle. In these iterations, the densities and change, 
but the interaction U, being a parameter of the Hamiltonian, does not. We have 
therefore expressed the integral in Eq.Q and the resulting transcendental equation 
for /3 exclusively in terms of U, in order to guarantee that the only slightly time- 
consuming steps of the calculation take place only once in each calculation, outside 
the selfconsistency cycle. This is the fifth additional condition on the spin- dependent 
parametrization, alluded to above. 

As a simple example of such KS calculations, which serves to illustrate all essential 
aspects, we consider open Hubbard chains, where the spatial inhomogeneity stems from 
the boundaries, which give rise to charge and spin-density oscillations in the bulk. 
Representative results for a chain with L = 100 sites are displayed in Fig. HJ where 
we compare the ground-state density profile (FigHJ-a) and spin-density profile (FigJU- 
b) obtained from DMRG to LSDA profiles obtained from using the fully numerical 
solution to the BA integral equations and from the present parametrization. The 
BALSDA/FVC and BALSDA/FN ground-state energies of the same system deviate 
from the DMRG energy by 0.01% and 0.64%, respectively. The local densities follow the 
same trend, and deviate from DMRG by 0.42% and 0.58%. For the local magnetization, 
the corresponding numbers are 2.20% for BALSDA/FVC and 5.86% for BALSDA/FN. 
Remarkably, for all three quantities the parametrized results are closer to the DMRG 
benchmark data than the numerically defined BALDA. This shows that the particular 
form chosen for the proposed parametrizations allows for considerable error cancellation. 
On 32 processors the DMRG calculation took approximately 17 hours, while the 
BALSDA calculations required approximately 40 seconds. Of course, for high-precison 
calculations, as well as for the calculation of quantities that are not easily extracted 
from densities and energies, DMRG is still essential. 

A more complex case is depicted in Fig. |5l which shows density and magnetization 
profiles for parabolically confined systems in a periodic chain. For two different values 
[k = 0.05 and k = 0.5) of the curvature of the confining potential (whose form is 
schematically indicated by the dashed (green) curve), the data points show the site- 
resolved particle density and spin density obtained by exact (Lanczcos) diagonalization, 
fully numerical BA-LSDA and our presently proposed parametrization. To within the 
accuracy that can be expected from a local-density approximation for this type of system 
(a few percent) the agreement between all three sets of calculations is excelent, for both 
charge and spin distributions. The amplitude of the magnetization-density oscillations is 
overestimated by both flavours of local approximations, which is consistent with previous 
observations for similar approximations and systems |20] . 

Next, we compare, in Fig. |6l the numerically exact ground-state energy obtained 
from Lanczos diagonalization of a small open Hubbard chain with L = 15 sites to the 
LSDA ground-state energies obtained from using the fully numerical solution to the 
BA integral equations and from using the present parametrization. Up to [/ ~ 4it the 
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Figure 4. (a) Density profile rii of an open cliain with L = 100 sites, = 30, = 20 
fermions, and U = it, obtained selfconsistently from fully numerical BALSDA, the 
present parametrization and DMRG calculations, (b) Local magnetization to,; of the 
same system. Percentage deviations, summed over all sites, of the LSDA densities from 
the DMRG ones are given in parenthesis. 




1 3 5 7 9 11 13 15 1 3 5 7 9 11 13 15 

site site 



Figure 5. Parabolically confined systems: (a) Density profiles Ui of periodic chains 
for two different values of curvature (A: = 0.05 and k — 0.5). L = 15 sites, Nf — 4, 
Ni = 3 fermions, and U — At, obtained selfconsistently from fully numerical BALSDA, 
the present parametrization and exact (Lanczos) calculations, (b) Local magnetization 
rrii of the same systems. 

BALSDA/FN data are almost identical to the exact data, which attests to the quality 
of the local approximation. For U larger than ~ 5t, the present parametrization is, 
once again, better than the conceptually superior fully numerical LSDA, due to error 
cancellation. The inset shows that for larger systems the behaviour is qualitatively the 
same. 

Finally, we point out that successful applications of our (then unpublished) spin- 
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Figure 6. Per-site ground-state energy of finite open chains with L = 15 sites, Nf = 4 
and Ni — 3 fermions, obtained selfconsistently from fully numerical BALSDA, the 
present parametrization and from exact diagonalization. Inset: BALSDA data for a 
larger chain with L = 200 sites, = 60 and = 40. 

dependent expression (|TT|) to the study of spin-polarized transport across a correlated 
nanoconstriction [11] and to the calculation of occupation probabilities of exotic 
superfiuids in spin-imbalanced systems [25] have aheady been reported. 

5. Summary 

In summary, we have constructed a simple and rehable parametrization for the ground- 
state energy of the homogeneous one- dimensional Hubbard model, for arbitrary filhngs, 
spin-polarizations and interactions. For the first time, a qualitatively and quantitatively 
correct description of the Mott gap is obtained from a simple density functional with 
a proper explicit derivative discontinuity. This parametrization can be used in its own 
right, for the homogeneous model, whenever simple expressions for the ground-state 
energy and for the resulting Mott gap are required. 

However, its main application is as input for local-density and local-spin-density 
approximations, which allow one to efficiently minimize the energy and extract energies, 
density profiles and related quantities for spatially inhomogeneous models. Since in KS 
calculations one never diagonalizes the interacting Hamiltonian, but only the auxiliary 
noninteracting one, systems with thousands of sites can be dealt with, even in the 
absence of any simplifying symmetry. 

This work was supported by Brazilian agencies CAPES, CNPq and FAPESP. The 
authors thank Dominik Horndlein for the DMRG data and Vivaldo Campo for his code 
for solving the Bethe-Ansatz integral equations and the numerically defined BA-L(S)DA. 
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